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Abstract.  The  Direct  Simulation  Monte  Carlo  (DSMC)  method  has  been  very  successful  for  the  study  of  many  problems 
in  rarefied  gas  dynamics  and  hypersonic  flow.  The  extension  to  applications  such  as  acoustics  will  provide  a  useful  tool 
for  capturing  all  physical  properties  of  interest  for  nonlinear  acoustic  problems,  such  as  dispersion,  attenuation,  harmonic 
generation  and  nonequilibrium  effects.  The  validity  of  DSMC  for  the  entire  range  of  Knudsen  numbers  ( Kn ),  where  Kn 
is  defined  as  the  mean  free  path  divided  by  the  wavelength,  allows  for  the  exploration  of  sound  propagation  at  low  Kn 
(low  frequency,  atmospheric  conditions)  as  well  as  sound  propagation  at  high  Kn  (high  frequency,  dilute  gases,  or  in 
microdevices).  This  paper  will  present  the  application  of  DSMC  to  nonlinear  acoustics  in  monatomic  and  diatomic  gases 
for  varying  values  of  Kn. 


INTRODUCTION 

For  high  amplitude  sound  propagation  at  atmospheric  conditions,  continuum  equations  such  as  the  nonlinear 
Euler  and  Navier-Stokes  equations  have  been  widely  used  to  understand  the  physical  properties  of  sound  due  to 
nonlinearity[l][2].  However,  an  overview  of  such  approaches  reveals  that  nonlinear  effects  are  often  treated 
separately,  and  then  superimposed  to  produce  a  final  result.  Instead  of  employing  this  technique,  a  method  such  as 
DSMC  can  model  all  nonlinear  and  viscous  effects  in  the  same  method.  Various  DSMC  acoustic  simulations  have 
been  compared  to  solutions  to  the  linear  and  nonlinear  Euler  equations  for  monatomic  gases  to  study  effects  such  as 
harmonic  generation,  absorption,  wave  steepening,  and  nonequilibrium  effects  in  a  previous  paper  [3]. 

The  Knudsen  number  {Kn)  is  defined  to  be  the  mean  free  path  divided  by  a  characteristic  length  scale  [4].  In  this 
study,  the  characteristic  length  scale  is  the  acoustic  wavelength  of  the  sound  wave.  Therefore,  the  Knudsen  number 
is  large  for  sound  propagation  in  very  dilute  gases,  in  micro-channels,  or  at  high  frequencies,  and  thus  requires  a 
particle  method  (or  kinetic  theory)  solution.  There  have  been  numerous  attempts  to  study  sound  propagation  for 
high  Kn  based  on  kinetic  theory.  Most  of  these  attempts  use  approximations  to  the  Boltzmann  collision  integral. 
For  example,  Wang  Chang  and  Uhlenbeck  [5]  studied  sound  propagation  using  the  method  of  eigenfunctions. 
Numerical  procedures  applied  to  the  eigenfunction  method  were  performed  by  Pekeris  et  al  [6].  In  Sirovich  and 
Thurber  [7],  sound  propagation  was  studied  using  a  kinetic  model  description  based  on  the  Bhatnagar  -  Gross  - 
Krook  (BGK)  equation  [8],  while  Marques  [9]  replaced  the  collision  integral  with  a  relaxation-time  term.  All 
attempts  report  that  the  speed  and  absorption  of  sound  depend  heavily  on  Kn.  Experimental  results  by  Greenspan 
[10],  Meyer  and  Sessler  [11],  and  Schotter  [12]  show  close  agreement  for  the  numerical  methods  above  for  Kn<  1, 
but  for  larger  Kn ,  the  results  are  poor. 

The  Direct  Simulation  Monte  Carlo  (DSMC)  [4,  13,  14]  method  is  a  particle  method  that  was  chosen  for  the 
current  study.  DSMC  describes  gas  flows  through  direct  physical  modeling  of  particle  motions  and  collisions. 
DSMC  was  originally  designed  for  aerospace  applications  in  dilute  gases,  and  thus  has  been  thoroughly  tested  for 
high  Kn  flows  (>0.2).  It  has  been  applied  to  hypersonic  and  rarefied  gas  flows,  and  has  in  fact  become  the  principle 
method  for  high  Kn  simulations  [15-18].  DSMC  has  also  been  successfully  applied  to  modeling  chemical  reactions 
[19-21],  simluating  blast-structure  interactions  at  low  Kn  [22],  and  to  study  sound  in  dilute  gases  (high  Kn)  [23]. 
However,  there  has  been  limited  work  in  using  DSMC  for  acoustics  and  low  Kn  problems  (continuum  regime). 

The  motivation  for  this  study  is  three  fold.  First  we  would  like  to  show  the  validity  of  using  DSMC  for  low  Kn 
applications.  Second,  we  would  like  a  numerical  method  to  include  all  viscous  and  nonlinear  acoustic  effects,  and 
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third,  we  would  like  a  method  to  accurately  model  sound  at  high  Kn.  A  parallel,  one  dimensional  implementation  of 
DSMC  for  monatomic  and  diatomic  gases  has  been  applied  and  compared  to  nonlinear  Euler  equation  solutions  as 
well  as  nonlinear  acoustic  theory  for  varying  values  of  Kn.  There  are  also  comparisons  to  experimental  data. 


DSMC  IMPLEMENTATION 

DSMC  is  based  on  the  kinetic  theory  of  gases,  where  representative  particles  are  followed  as  they  move  and 
collide  with  other  particles  [4].  The  movement  of  particles  is  determined  by  their  velocities,  while  the  collisions  are 
determined  statistically,  but  are  required  to  satisfy  mass,  momentum  and  energy  conservation.  DSMC  is  usually 
limited  to  binary  collisions.  In  DSMC  the  particle  positions  and  velocities  are  initialized  randomly  and  the  boundary 
conditions  ultimately  determine  the  final  solution,  as  in  continuum  methods. 

DSMC  simulations  become  more  and  more  exact  as  the  time  step  and  the  cell  size  tend  to  zero,  and  has  in  fact 
been  proven  to  converge  to  the  Boltzmann  equation  in  the  limit  of  small  time  step  and  small  cell  size  [24].  It  can  be 
shown  that  the  Boltzmann  equation  will  reduce  to  the  Navier-Stokes  equations  for  Kn<  0.1  and  the  Euler  equations 
for  Kn  ->  0.  Therefore,  at  low  Kn,  this  means  that  DSMC  simulations  and  Euler  equation  solutions  should  give  the 
same  results.  Also,  the  Boltzmann  equation  includes  all  physical  properties  of  interest  for  sound  propagation. 
Without  modification  to  the  method,  properties  such  as  viscous  dissipation,  nonlinear  effects,  dispersion, 
nonequilibrium  effects,  entropy  waves,  and  more  are  available.  In  fact,  DSMC  has  been  able  to  model  more 
physical  phenomena  than  is  usually  contained  within  the  Boltzmann  equation,  such  as  chemical  reactions,  multiple 
species,  and  diatomic  particles. 

In  this  study,  acoustic  waves  were  simulated  by  creating  a  pressure  distribution  in  the  middle  of  the  domain, 
which  then  propagates  in  both  directions.  The  shape  of  this  pressure  distribution  is  dependent  on  the  acoustic  source 
of  interest.  Here,  Gaussian  and  sinusoidal  pressure  distributions  of  varying  frequencies  were  considered.  The  flow 
velocity  is  initially  set  to  0.0,  and  both  monatomic  and  diatomic  hard  sphere  molecules  with  binary  elastic  collisions 
were  simulated.  The  cell  sizes  were  always  less  than  !4  a  mean  free  path,  and  care  was  taken  in  the  choice  of  At  to 
ensure  accuracy,  which  was  computed  by  using  a  Courant  or  von  Neumann  number  of  0.1  in  most  cases.  At  must 

also  be  less  the  mean  collision  time.  This  corresponds  to  time  step  sizes  in  the  range:  10  12  <  At  <  10  11 . 


COMPUTATION 

An  Object  Oriented  Programming  (OOP)  [25,  26]  approach,  written  in  C++,  is  used  in  the  development  of  this 
DSMC  study.  Such  an  implementation  allows  the  code  to  be  more  readable,  manageable,  reusable,  and  extendable. 
The  program  is  also  more  efficient  by  using  such  C++  characteristics  such  as  memory  handling  and  pointers.  Using 
the  object  oriented  C++  concepts  of  class  structure ,  the  code  is  easy  to  modify  and  understand.  An  object  oriented 
technique  is  a  very  natural  approach  for  this  problem  because  the  particles  and  cells  are  physical  objects.  In 
addition,  OOP  allows  encapsulation,  the  grouping  together  of  data  and  methods  (i.e.  functions).  Using  objects,  the 
direct  simulation  Monte  Carlo  method  flows  in  a  logical  sequence  and  is  easy  to  understand.  The  responsibility  of 
the  detailed  calculations  is  contained  in  the  individual  object. 

Since  the  code  runs  independent  simulations  and  then  performs  ensemble  averages  to  reduce  scatter,  it  lends 
itself  very  nicely  to  a  parallel  computing  algorithm.  Instead  of  having  one  processor  do  all  the  independent  runs,  it 
is  easy  to  have  different  computers  simultaneously  run  different  ensembles,  and  then  average  them. 

To  achieve  this  simple  parallel  processing  algorithm,  a  master/slave  routine  is  constructed.  The  master’s  main 
job  is  to  schedule  the  processors  for  work,  organize  collected  data,  and  do  the  final  computing  of  needed  flow 
properties.  The  slaves  compute  the  ensemble  runs  and  send  sampled  data  back  to  the  master  for  final  processing. 
The  programming  here  is  made  possible  by  using  the  Message  Passing  Interface  (MPI)  [27]  to  send  messages  back 
and  forth  between  processors.  Also,  the  programming  for  this  type  of  algorithm  is  quite  simple  because  there  are 
very  few  MPI  calls. 

The  code  was  run  on  a  160-processor  Beowulf  cluster  running  Linux.  It  has  a  Dolphin  network  for  fast 
communication  and  2.8  GHz  Athlon  processors.  It  has  a  peak  speed  of  0.4  teraflops  (floating  point  operations  per 
second).  Despite  the  efficiency  of  the  parallel  algorithm,  the  CPU  time  for  these  DSMC  programs  is  still  fairly  large 
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ACOUSTICS  THEORY  AND  RESULTS 


For  a  sound  wave  with  high  amplitude,  the  linear  wave  equation  cannot  accurately  simulate  wave  propagation. 
Nonlinear  acoustic  effects  will  therefore  become  more  evident  as  the  amplitude  becomes  large.  Some  of  the  major 
nonlinear  acoustic  effects  are  wave  steepening,  shock  formation,  harmonic  generation,  and  dissipation  [1,  2].  The 
speed  of  sound  for  a  perfect  gas  is  given  to  be: 

c  =  pRT  (i) 


which  will  be  different  for  monatomic  and  diatomic  gases.  Because  of  the  large  acoustic  pressures  in  high 
amplitude  sound,  pressure  peaks  travel  faster  than  pressure  troughs.  An  implication  of  this  is  that  the  wave  steepens 
with  increasing  time  and  propagation  distance. 

Steepening  continues  until  the  sound  wave  begins  to  shock.  This  happens  at  a  distance  X  which  is  defined  to  be: 

2 


X  = 


A)  co 

j3kP0 


(2) 


Where  pQ  is  the  ambient  density,  C0  is  the  ambient  sound  speed,  k  is  the  acoustic  wave  number,  PQ  is  the  initial 
acoustic  pressure  amplitude,  and  j3  is  the  coefficient  of  nonlinearity,  which  has  the  value  of  {j  +  l)/2  .  This  gives 

/?  =  4/3  in  a  monatomic  gas  and  /3  =  1 .2  in  a  diatomic  gas.  For  low  frequency,  low  amplitude  sound  waves,  this 
shock  distance  is  typically  very  large,  where  the  signal  will  dissipate  before  the  wave  is  given  a  chance  to  shock. 
However,  for  high  frequency,  high  amplitude  sound  waves,  a  shock  wave  can  develop  within  a  few  wavelengths. 
For  sound  in  diatomic  gases,  Equation  (2)  implies  that  the  shock  distance  will  be  greater  than  the  same  frequency  in 
a  monatomic  gas. 

DSMC  simulations  were  implemented  in  monatomic  and  diatomic  gases  to  investigate  nonlinear  and  viscous 
effects.  The  diatomic  routine  uses  the  phenomenological  Larsen-Borgnakke  method  [4,  28]  to  compute  internal 
energy.  Results  for  a  DSMC  simulation  of  a  Gaussian  pulse  in  Oxygen  are  shown  in  Figure  1  after  it  has  propagated 
some  distance.  Notice  that  the  Gaussian  pulse  has  steepened  significantly  and  the  signal  has  become  distorted  due  to 
nonlinearity.  Figure  2  also  shows  the  differences  between  monatomic  simulation  in  Argon  and  a  diatomic 
simulation  in  Oxygen.  It  appears  that  the  diatomic  case  has  not  steepened  as  much  as  the  monatomic  case  which  we 
would  expect  since  the  shock  distance  is  further  for  Oxygen. 


FIGURE  1.  Comparisons  of  DSMC  results  in  Oxygen  against  the  nonlinear  Euler  equations. 
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FIGURE  2.  Comparisons  of  DSMC  results  in  monatomic  and  diatomic  gases 

The  DSMC  results  were  compared  to  a  finite  difference  solution  of  the  nonlinear  Euler  equations  that  uses  a  time 
accurate  McCormack  scheme  and  added  artificial  viscosity  to  reduce  numerical  oscillations.  The  Gaussian  pulse  is 
100  cells  wide  with  each  cell  being  lA  of  a  mean-free-path  wide.  The  mean-free-path  in  these  conditions  is 
approximately  50  nanometers.  The  domain  is  about  1000  cells  long.  Each  cell  has  about  150  particles.  There  were 
4000  time  steps,  each  of  size  7.5  picoseconds.  The  acoustic  pressure  amplitude  of  the  Gaussian  source  is  equivalent 
to  190  dB,  with  acoustic  flow  velocity  equal  to  zero.  Here,  Kn  =  0.02.  It  is  clear  that  the  comparison  between 
DSMC  and  nonlinear  Euler  solution  in  this  case  gives  excellent  results. 

Due  to  wave  steepening,  energy  from  the  fundamental  frequency  gets  shifted  to  the  higher  harmonics,  allowing 
the  harmonics  amplitudes  to  grow  and  thus  distorting  the  signal.  The  harmonic  component  amplitude  can  be 
determined  analytically  for  a  sine  wave  in  a  lossless  medium  by  the  Fubini  solution  [1,2]  and  the  Blackstock 
Bridging  function  [1].  DSMC  results  for  a  monatomic  gas  are  compare  to  the  Blackstock  bridging  function  in 
Reference  [3].  The  Fourier  component  amplitude  of  the  first  4  harmonics  of  a  30  MHz  DSMC  sine  wave  were 
compared  to  the  analytic  Blackstock  Bridging  function.  The  agreement  was  excellent. 

Rarefied  Gas  Acoustics 

Experimental  and  other  kinetic  theory  results  suggest  that  the  absorption  of  sound  for  high  Kn  differs 
significantly  from  that  for  low  Kn  [5-7,  9,  29].  The  scaled  absorption  coefficient,  a/k 0  ,  where  a  is  the  absorption 

coefficient  and  £0is  the  acoustic  wave  number,  for  Kn  ranging  from  0.04  to  35  is  shown  in  Figure  3  for  both 

monatomic  and  diatomic  cases.  DSMC  simulation  results  are  compared  to  a  classical  Navier-Stokes  prediction, 
experimental  work  by  Greenspan  [10],  and  previous  kinetic  theory  results  by  Hadjiconstantinou  and  Garcia  [23] 
whose  work  was  based  on  standing  waves  in  waveguides.  The  variation  in  Kn  was  obtained  by  keeping  the  cells  per 

wavelength  constant  at  100,  but  varying  the  cell  size  from  1/2  of  a  mean  free  path  to  3/1 000  of  a  mean  free  path. 

The  absorption  coefficient  results  were  obtained  by  finding  the  maximum  pressure  at  each  time  step  and  calculating 
a  line  of  best  fit  on  a  logarithmic  scale. 
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FIGURE  3.  Scaled  absorption  coefficient  for  varying  Kn  ranging  from  0.04  to  35 


For  high  Kn  simulations,  the  wavelength  of  the  signal  is  significantly  shorter  than  a  mean  free  path. 
Experimental  results  for  Kn  >  5.5  are  thus  difficult  to  obtain  because  of  scattering  and  high  absorption.  The 
simulations  presented  here  have  successfully  captured  the  physics  for  such  cases  by  decreasing  our  cell  size  to 

3/1000  a  mean  free  path.  The  calculations  confirm  that  the  Navier-Stokes  prediction  for  the  absorption  coefficient 
is  only  correct  for  low  Kn,  and  a  kinetic  theory  approach  must  be  used  for  high  Kn  calculations.  The  absorption 
coefficient  results  are  in  good  agreement  with  experimental  results  and  previous  work  on  standing  waves. 

High  Kn  acoustics  can  also  involve  nonequilibrium  effects,  which  is  another  reason  a  kinetic  theory  approach  is 
needed  to  study  sound  at  high  Kn.  Figures  4  and  5  show  the  velocity  distribution  functions  for  Kn  =  0.02  and  5  in 
Oxygen.  Both  are  compared  to  the  Maxwellian  distribution,  and  both  are  taken  at  the  pressure  peaks  of  the  sound 
wave.  For  Kn  =  0.02,  the  distribution  fits  the  Maxwellian  very  well,  as  we  would  expect  this  system  to  be  at 
equilibrium.  However,  for  Kn  =  5,  there  is  a  shift  from  equilibrium,  as  noticed  in  Figure  5. 


FIGURE  4.  Velocity  distribution  curves  at  Kn  =  0.02  FIGURE  5.  Velocity  distribution  curves  dXKn  =  5 


557 


CONCLUSIONS 


The  current  study  uses  DSMC  to  study  nonlinear  acoustic  propagation.  This  model  was  implemented  using  object 
oriented  programming  in  C++  which  makes  the  entire  code  easier  to  read,  maintain,  and  expand.  The  code  also  uses 
a  parallel  processing  algorithm  to  further  improve  the  efficiency  of  the  model.  Using  DSMC  for  the  study  of 
nonlinear  acoustics  was  validated  by  comparisons  to  the  Euler  equations,  the  Navier-Stokes  equations,  and 
experiment.  Sound  propagation  for  high  Kn  is  discussed  where  the  absorption  of  sound  is  higher  than  a  Navier- 
Stokes  or  Euler  prediction.  Good  agreement  is  found  between  experiment  and  previous  work.  Results  indicate  that 
DSMC  is  a  suitable  technique  for  modeling  sound  propagation  for  all  values  of  Kn  for  both  monatomic  and  diatomic 
gases. 
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